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We derive a first order differential equation for the decoherence of an orbital angular momentum 
entangled biphoton state propagating through a turbulent atmosphere. The derivation is based 
on the distortion that orbital angular momentum states experience due to propagation through a 
thin sheet of turbulent atmosphere. This distortion is treated as an infinitesimal transformation 
leading to a first order differential equation, which we call an infinitesimal propagation equation. 
The equation is applied to a simple qubit case to show how the entanglement decays. 
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I. INTRODUCTION 

The orbital angular momentum (OAM) states of a pho- 
ton is an attractive basis for quantum information pro- 
cessing and communication, because it allows a higher 
dimensional representation of quantum information [lj- 
3]. One application where such a higher dimensional 
representation can have benefits in terms of information 
capacity is in free-space quantum communication. Un- 
fortunately, the turbulence in the atmosphere causes the 
decoherence of entanglement of OAM states. The refrac- 
tive index fluctuations in the turbulent atmosphere pro- 
duce random phase modulations that distort the OAM 
states. This process is similar to the distortion suffered 
by OAM modes in classical optical beams due to scintil- 
lation in turbulence 



In this paper we derive a theoretical framework with 
which one can investigate the decoherence of OAM en- 
tanglement in atmospheric turbulence. Some previous in- 
vestigations of this process are based on the work 
of Patterson who assumed that one can model the 
turbulent medium with a single phase screen and param- 
eterized the turbulence with one parameter, the Fried 
parameter [12]. As a result the effects of the propaga- 
tion, such as beam spreading, intensity scintillation and 
multiple scattering among different modes, are not incor- 
porated into the model. Neither is this framework able 
to consider the individual effects of the different scale 
parameters that are combined into the Fried parameter. 

Another approach to investigating the decoherence 
of OAM entanglement in atmospheric turbulence is to 
model the turbulent medium as an absorbing and scat- 
tering medium [l3j. However, this analysis assumes that 
the turbulence is weak. 

The approach that is used in this paper is based on 
the incremental effect of the turbulent atmosphere on the 
quantum state (density operator) , expressed as a first or- 



der differential equation, which can then be solved to ob- 
tain the evolution over the entire propagation distance. 
This approach is reminiscent of a master equation ap- 
proach, however, here, the derivative is taken with re- 
spect to propagation distance instead of time. In this 
framework one can incorporate any model for the tur- 
bulence — i.e. any power spectral density such as Kol- 
mogorov, Tartaskii or von Karman By implication 
one can investigate the effect of any of the scale param- 
eters associated with the whole process and not only the 
Fried parameter. One can apply this framework to cases 
with arbitrary strong turbulence or strong scintillation. 

For the purpose of the derivation, it is assumed that 
the light is monochromatic, that the beam propagates 
paraxially and that it has a uniform polarization to allow 
scalar propagation. The final result can be generalized to 
polychromatic vector fields. The paraxial approximation 
would always be valid in practical applications of this 
framework. 

The paper is organized as follows. First we discuss how 
to treat a photon field that propagates through a lin- 
ear spatially inhomogeneous medium in Section [TTJ Then 
we discuss the infinitesimal transformation of the mo- 
mentum space wave functions of such photon fields in 
Section IIII1 In Section IIVI we show how we define the 
density operator in terms of the OAM and momentum 
bases. The derivation of the infinitesimal propagation 
equation (IPE) is shown in Section [Vj with the aid of 
the ensemble averaging discussed in Appendix [Bj Then 
in Section IVII we show how one can treat the different 
integrals that appear in the IPE, using the generating 
function of OAM model discussion in Appendix \K\ In 
Section I VIII a simple example is considered. We discuss 
the results in Section IVIIII and end with a summary in 
Section IS 



II. FIELD QUANTIZATION 
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A dynamical system is generally understood to be one 
evolving in time. For a quantum system this evolution 
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is unitary, being described by a unitary operator that is 
given by an exponential operator with the time integrated 
Hamiltonian as its argument. Within this context the 
process of decoherence of OAM entanglement in a turbu- 
lent atmosphere is not described by a dynamical system. 
Although the turbulent medium is inhomogeneous it is 
still assumed to be linear (l7| . which implies that one 
can consider each temporal frequency separately, leading 
to a monochromatic assumption. The temporal behav- 
ior of the field is therefore simply given by exp(zu;£). A 
straight forward Hamiltonian approach for this problem 
is therefore inappropriate [l8| . Instead of having a three- 
dimensional field that evolves in time, our approach in 
this paper is to ignore the temporal behavior and instead 
consider a two-dimensional field that changes as it moves 
along the third spatial dimension. 

Due to the monochromatic assumption, which fixes the 
wavelength of the light A, the ^-component of the propa- 
gation vector can be expressed as a function of the other 
two components 

/47T 2 \ 1 / 2 

k z (k x , k y ) = I k x — k y j . (1) 

As a result the momentum basis under the monochro- 
matic approximation becomes a two-dimensional basis 
|K), where the propagation vector K = k x x + k y y 
represents the two-dimensional projection of the three- 
dimensional propagation vector. It is necessary to in- 
clude the evanescent momentum states (those for which 
k z is imaginary) to ensure that the basis is complete. 
However, under the paraxial approximation the evanes- 
cent part of the field is negligible. We assume that the z- 
direction is the general direction of propagation. There- 
fore, the two-dimensional momentum basis describes two- 
dimensional fields on the transverse plane, perpendicular 
to the propagation direction. In the absence of turbu- 
lence one can recover the full three-dimensional field by 
using Fresnel diffraction theory to propagate the two- 
dimensional field along z. In the presence of turbulence 
the two-dimensional field transforms in a more compli- 
cated way from one plane to another. We formulate the 
decoherence process in terms of this transformation. 

An arbitrary pure state of a monochromatic, uniformly 
polarized single photon can now be expressed on the 
transverse plane in terms of the two-dimensional momen- 
tum basis 

W) = J G(K)|K) 0, (2) 

where the coefficient function G(K) is the momentum 
space wave function (K|?/>). The inverse two-dimensional 
Fourier transform of the momentum space wave func- 
tion gives the position space wave function g(x, y) on the 
transverse plane at z = 0. 

Although the OAM basis (which is equivalent to the 
two-dimensional momentum basis) is the preferred basis 
for an analysis where the photons are entangled in term 



of OAM, the turbulence is modeled in terms of a power 
spectral density that is defined in terms of a momentum 
basis. As a result one is forced to have a mixture of 
both bases in this analysis. We therefore start with the 
momentum basis and eventually express the OAM basis 
in terms of it. 



III. EQUATION OF MOTION AND 
INFINITESIMAL TRANSFORMATION 

For an interaction-free system the quantum wave func- 
tion obeys the same equations of motion as the classical 
field. Therefore, the equation of motion for the position 
space wave function can be derived from Maxwell's equa- 
tions. In a source- free region this gives the Hclmholtz 
equation 

V 2 £(x) + n 2 k 2 E(x) = 0, (3) 

where E(x) is the scalar part of the electric field (as- 
suming the polarization is uniform and can be ignored), 
k is the wave number, n is the refractive index and 
x = xi + yy + zz. The inhomogeneous medium is repre- 
sented by a spatially varying index of refraction 

n=l + 5n(x). (4) 

This variation is very small (Sn <C 1), which implies that 
one can approximate the Helmholtz equation as 

V 2 E(x) + k 2 E(x) + 2Sn(x)k 2 E(x) = 0. (5) 

Furthermore, we assume that the beam is paraxial and 
propagates in the z-direction. So we define 

E(x) = ,g(x) exp(-ifo), (6) 

which then leads to the paraxial wave equation with the 
extra inhomogeneous medium term 

V|g(x) - i2kd z g(x) + 2Sri(x)k 2 g(x) = 0, (7) 

where Vj is the transverse part of the gradient operator. 

Now we replace <?(x) with its two-dimensional inverse 
Fourier transform, 

5 (x) = J G(K, z) cxpHK • x) (8) 

which contains the angular spectrum of the optical field 
G(K, z). The latter also represents the momentum space 
wave function for the purpose of the quantum analysis. 
Then we obtain 

d z G{K,z) = ^\K\ 2 G(K,z)-ikN{K,z)*G{K,z), (9) 

where * indicates convolution and N(K, z) is the two- 
dimensional Fourier transform of <5n(x). This equation 
represents the infinitesimal transformation of the mo- 
mentum space wave function during propagation. It 
forms the basis of the derivation of the IPE. 
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IV. DENSITY OPERATOR IN THE OAM BASIS 

For the purpose of this derivation we first consider a 
single photon and then generalize the result for the case 
of two photons. The density operator of an arbitrary 
single photon state can be expressed in the OAM base 
by 



m) p mn {n\ 



(10) 



where we use a single index to denote both the indices 
of the OAM states, m = {l,r}. (See Appendix lAl for 
a discussion of the OAM modes.) Unless stated other- 
wise, each OAM index used in the subsequent deriva- 
tion always represents both the indices associated with 
a particular OAM mode. Each of these OAM states can 
be expanded in terms of the two-dimensional momentum 
basis 



G m (K)\K) 



d 2 K 
Air 2 ' 



(11) 



leading to the following expression for the density oper- 
ator in terms of the two-dimensional momentum basis 



Pi?) 



/ |Ki) G m (Kx,z)p mi7l 
d 2 K 1 d 2 A 2 



xg;(k 2 ,z)(k 2 | 



47T 2 47T 2 



(12) 



where the dependence on z is shown explicitly to make 
it apparent that this expression for p is only valid on 
a transverse plane for a specific value of z. Here, and 
also later, we use only one integral sign to represent sev- 
eral A-space integrals. By evaluating the summations in 
Eq. ([i"2j) . one obtains the definition for the density oper- 
ator in terms of the two-dimensional momentum basis, 



where 

p(Ki,K 2 ,z 

which implies that 



J |Ki)p(Ki,K 2 ,z)(K 2 



d 2 Ai d 2 A 2 
4tt 2 4tt 2 



(13) 



m.n 



G m (K 1 ,z)p min G* n (K 2 ,z), (14) 



G* m (K 1 ,z) j0 (K 1 ,K 2 ,z)G„(K 2 ,z) 



47T 2 47T 2 

(15) 

Since the two-dimensional momentum basis and the 
OAM basis are completely equivalent, the definitions in 
Eqs. (fT0| and (fT3|) are also completely equivalent and 
Eqs. (I14[) and (|15[) indicate how one can transform from 
one to the other. 

For two photons the density operator in Eq. (|12|) can 
be generalized to become 



p(z) 



J |K 1 )|K 3 )G m (K 1 ,z)G p (K. 



3, Z) 



Xp v 



„ p<q G;(K 2 ,z)G*(K 4 ,z)(K 2 |(K 4 | 



d 2 Ki d 2 K 2 d 2 K 3 d 2 K 4 

47T 2 47T 2 47T 2 47T 2 



(16) 



To obtain the full three-dimensional expression for the 
density operator in free-space (without turbulence) one 
can use Fresnel diffraction theory to determine the ex- 
pression at any other value of z. In the presence of tur- 
bulence the expression for the density operator is only 
valid on a specific transverse plane, and it needs to be 
transformed according to Eq. ^) from plane to plane. 

Note that inside the expression the z-dependence is 
carried by the momentum space wave functions and not 
by the density matrix elements. This is because the 
transformation of the density operator during propaga- 
tion over an infinitesimal distance through turbulence is 
caused by the distortion of the momentum space wave 
functions. After such an infinitesimal propagation these 
momentum space wave functions no longer represent the 
Fourier transforms of the original modes. One needs 
to re-expand these distorted wave functions in terms of 
the momentum space wave functions of the OAM modes 
and incorporate the expansion coefficients in the den- 
sity matrix elements. Thereby one can transfer the z- 
dependence to the density matrix elements. 
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FIG. 1: Diagram showing how the OAM entangled biphoton 
propagates through turbulent media toward the two detec- 
tors. 



For notational convenience we represent the product of 
momentum space wave functions that appear in Eq. (|16|) , 
as a single function, 

G ro (K 1) z)G p (K 3) z)G* n (K 2 , z)G* q (K 4 , z) 
= F rWi9 (K 1 ,K 2 ,K 3 ,K 4 ,z). (17) 



d 2 K\ d 2 A 2 The expression for p is now given by 



p(z) 



J2 J |Ki)|K 3 ) Pn 



J m,n,p,q 

m,n,p,q ' 

Xim,tl,J),g(Kl) K 2 , K3, K4, z) 

^3 



x(K 2 |(K 4 ; 



d 2 Ai d 2 A 2 d 2 A 3 d 2 A 4 

47T 2 47T 2 47T 2 47T 2 



(18) 



Since F m ^ n ^ p ^ q carries the only z-dependence in the ex- 
pression for the density operator, we focus on how it 
transforms during infinitesimal propagation. At the end 
we apply the transformation to the expression for p. 



DERIVATION OF THE IPE A. Infinitesimal transformation 



We are ready to consider the steps in the derivation 
of the IPE. For this purpose we consider the scenario 
where a pair of photons (biphoton), which are entangled 
in terms of the OAM basis, both propagate through a 
turbulent atmosphere, as shown in Fig. [TJ It is assumed 
that the turbulences seen by the respective photons are 
mutually uncorrelated. 



The transformation that is caused by an infinitesimal 
propagation is obtained by substituting Eq. © into the 
infinitesimally propagated version of F mn ^ vq , which is 
obtained by setting z = z + dz and then expanding the 
result to subleading order in dz. The resulting expression 
is turned into a differential equation by taking the limit 
dz ->■ 0, 



(K ll K 2 ,K 3 ,K4^) 



— (|Ki| 2 - |K 2 | 2 + |K 3 | 2 - |K 4 | 2 )F m ,„ iM (K 1 ,K 2 ,K 3 ,K4,z) 
-ik J N 1 (K 1 -K,z)F„ tq (K,K 2 ,K 3 ,K A ,z) 

-J 7Vr(K 2 -K,z)F TO ^(K 1 ,K,K 3 ,K 4 ,z) 

d 2 ^ 



+ J N 2 (K 3 -K,z)F m!n>Piq (K 1 ,K 2 ,K,K i ,z) ^ 
- J iV*(K 4 - K, 0)F ro , n , M (K 1) K 2 ,K 3! K, z) 



d 2 K 
4^2" 



(19) 



r 



This differential equation still contains the iV n 's, which 
are random functions. One needs to perform an ensemble 
averaging to remove the randomness. However, if one 
would perform such an ensemble averaging on Eq. (|19[) 
all the terms containing N n 's would fall away, because 
(N n ) — 0. So we proceed as follows. Integrating Eq. (fl9l) . 

in terms of previous versions of itself. 



one obtains F, 



m,n,p,q 



Next, one substitutes the resulting integral equation back 
into itself repeatedly to obtain an infinite series (Dyson 



expansion). The ensemble average of the result would 
remove many terms, because 

(N n ) = (JV„JVW) = for m ^ n. (20) 

The resulting integral equation, up to second order in 
N n , showing only the terms that survives the ensemble 
averaging, is given by 



^m,n,p,g(Kl, K 2 , K3, K4, z) — F mt n iPiq (Ki, K 2 , K3, K4, Zo) 

+ (z - z ) A(| Kl | 2 - |K 2 | 2 + |K 3 | 2 - |K 4 | 2 )F m: „, p , 9 (K 1 , K 2 , K 3 , K 4 , z Q ) 

-k 2 [ r J (JVx(Ki -K,z 1 )7V 1 (K-K ,z 2 ))F TO , 9 (K ,K 2 ,K 3 ,K 4 ,zo) 

-<JV x (Ki - K, Zl )N*(K 2 - K , z 2 ))F„ hn ^ q (K, K , K 3 , K 4 , z ) 
-(N*(K 2 - K, *i)JVi(Ki - K , z 2 ))F rmq (K , K, K 3 , K 4 , z Q ) 
+ (N*(K 2 - K, Zl )N*(K - Ko, 2 2 ))F m>n , p , g (K 1) K , K 3 , K 4 , z Q ) 
+ (N 2 (K 3 - K, Zl )N 2 (K - K , z 2 ))F m , n , p ,,(K 1 , K 2 , K , K 4 , z Q ) 
-(N 2 (K 3 - K, zi)7V*(K 4 - K , z 2 ))F m>ntP<q (K 1 , K 2) K, K , z ) 
-(iV 2 *(K 4 - K, z!)iV 2 (K 3 - K , z 2 ))F m ^ q {^¥. 2 , K , K, z ) 

+ (N 2 *(K 4 -K,z 1 )N 2 *(K-K ,z 2 ))F m ^ q (K 1 ,K 2 ,K 3 ,K Q ,z Q ) ^ dz 2 d?i 

(21) 
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Note that the only z\- and ^-dependences appear in the N n 's. In Appendix [B] it is shown that 

f r (N n (K 1> z 2 )N*(K 2> z 1 )) dz 2 dz! = 2vr 2 dz<5(K 1 - K^^KO, (22) 

J za J Zn 



where we set z — zq = dz. We now use Eq. ([22]) to simplify Eq. (|2Tj) and then take the limit dz — > 0, to turn it into a 
differential equation again, 

0,J7 TO)niP , 9 (K 1 ,K 2 ,K3,K 4j aO = ^(|Ki| 2 - |K 2 | 2 + |K 3 | 2 - |K 4 | 2 )F mi „ iP)g (K 1) K 2l K 3 ,K4^o) 

d 2 X 



4tt 2 

d 2 ^ 



-2fc 2 F min , p , 9 (K 1 ,K 2 ,K 3 ,K4,zo) / * X (K) 
+k 2 J $ 1 (K)F myn ^ P:q (K 1 K, K 2 K, K 3 , K 4 , z ) 
+fc 2 y" $ 1 (K) J F m ,„^(K 1 ,K 2 ,K 3 -K,K4-K,zo) (23) 



4tt 2 

2 , 



If one substitutes Eq. (|23|) into the z-derivative of Eq. ([18]) one would obtain a first order differential equation for the 
density operator. However, we are interested in the transformation of the individual density matrix elements. 

B. Extraction of matrix elements 

To express the transformation of the density operator due to the infinitesimal propagation through a turbulent 
atmosphere in terms of the density matrix elements, one needs to extract the matrix elements from the density 
operator, using the trace 

d z pu,v,r, s (z) = trace {d z p(z)\v)\s)(u\(r\} , (24) 
where the operator that selects a particular matrix element in the OAM basis is given by 

\v)\s)(u\(r\ = J |K 8 )|Ka) F; w (K 5 , K 6 , K 7 , K 8 , z)(K 7 | (K 5 | (25) 

Substituting Eqs. (JUJ and flU) into Eq. ([24]) we obtain 

d z p u ,v,r,s(z) = traced ^ / |Ki)|K 3 ) Pm,n,p, g flg-Fm,n,p,g(Ki, K 2 , K 3 , K 4 , z)(K 2 \ (K 4 | "^F^F"!^ 

^ m,n,p,q 

d 2 K 5 d 2 K e d 2 K 7 d 2 K 8 ' 



x||k 8 )|k 6 ) f; w (k 5 ,k 6 ,k 7 ,k 8 ,z)(k 7 |(k 5| m 4 _ 2 , ;2 4/ J 

Pm 7 n,p,q J dzFm,n,p,q (Ki , K2 , K3 , K4 , z) 



m,n,p,q 



P * /"w \r \c \s ■d 2 K 1 d 2 K 2 d 2 K 3 d 2 K 4 
xF W| ,(K 1 ,K 2 ,K 3 ,K 4 ,z) ^ ^ ^ 2 ^ , (26) 

where the last expression is obtained because of the orthogonality of the momentum basis 

(K!|K 2 ) =4tt 2 ( 5 2 (K 1 -K 2 ), (27) 
as applied for the separate photons. Substituting Eq. ([2"B"[) into Eq. ([2l)]) . we obtain 



i 

2k { 



$zPu,v,r.s (^) — ^ ^ Pm.n,p,q i 
m,n,p,q 

-2fc 2 F m ,„ :Pi9 (K!, K 2 , K 3 , K 4 , z ) / *i(K) 



K 4 | 2 - |K 2 | 2 + |K 3 | 2 - |K 4 | 2 )F m ,„ iM (K 1 ,K 2 ,K 3 ,K 4 ,z ) 



d 2 K 

d 2 K 



+k 2 J $ 1 (K)F m>n , p , g (K 1 -K,K 2 -K,K3,K 4 ^o) ^ 



6 



d 2 K~ 



$ 1 (K)f , m , njM (Ki, K 2 , K 3 - K, K 4 - K, z ) ^ 

- J?* (XC K *T *T \ d2Rl d2R2 d2R3 d2Ri 

•\ w (Ki,K2,K3,K4,z) 4tt 2 4tt 2 4tt 2 4tt 2 ' 



(28) 



C. Final expression 



VI. SOLVING THE INTEGRALS 



Due to the orthogonality of the momentum space wave 
functions of the OAM basis, 



one can simplify Eq. (|28|) . The resulting first order dif- 
ferential equation, which represents the IPE, is given by 

&zPu,v,r,s(z) — $m,u\Z)Pm,v,r,s ^v,n\%) Pu,n,r,8 

iim,fi,u,i)(^)Pm,ti,r,s Lp^q^ TjS (z*)Pu,v,p,q 
— 2LTPu,v,r,s, (30) 

where repeated indices imply summation and the follow- 
ing definitions were made: 

s* tV (*) = 4/ i k i 2g *( k ^) g ;( k ^) Sr> ( 31 ) 



47T 2 



l t = r 



$i(K) 



d^K 
4tt 2 ' 



and 



W.vW = fc2 / $i(K)W m , u (K, z)VK„* „(K, z) 



(32) 



d 2 A- 



4tt 2 ' 
(33) 
with 

W X)J ,(K,z) = J G x (K 1 ,z)G* y (K 1 --K,z) (34) 

The first four terms of Eq. ([50)1 are the non-dissipative 
terms, representing the free-space propagation process. 
The last three terms of Eq. (f5T7|) are the dissipative terms, 
representing the scattering of OAM modes into other 
OAM modes due to the turbulence. One can separate 
the dissipative terms into separate pairs for the two pho- 
tons, each having its and Lt terms. 

The IPE in Eq. §U\), together with Eqs. (j3"T) -([M l) . 
is the main result of this paper. In general, Eq. (|30|) 
represents an infinite set of coupled first-order differen- 
tial equations. Even if the initial state contains only a 
few lower order modes, the turbulence will cause these 
modes to be coupled into all other modes. Subsequently, 
the other modes will couple back into the original modes. 
Truncating the set of equations, one inevitably excludes 
part of the coupling among all the different modes. How- 
ever, this coupling should become progressively smaller 
for higher order modes. Hence, one may be able to trun- 
cate the set at some point while retaining the dominant 
inter-modal coupling. 



The momentum space generating function, discussed 
in Appendix [XJ and given by 



J G m (K, z)G* n (K, z)^= 6 m , n , (29) F{G} 



1 



exp 



iir(a + ib)p iir(a — ib)q 



1 



1 



7T 2 (a 2 + b 2 )VL{t,w) 



- w 



(35) 



is now used to evaluate the integrals for Eqs. (|3T|) - ([M)) . 

The formalism allows one to use any power spectral 
density $o(k) for the turbulence. Here we neglect the 
effect of the inner scales, and use the von Karman power 
spectral density fl2j . 



$ (k) 



0.033C 2 



(|k| 2 + K 2 )H/6 



*l(*0, 



(36) 



where C 2 is the refractive index structure constant for 
the turbulence and n is inversely proportional to the 
outer scale of the turbulence. The outer scale helps to 
regularize the integrals, but in the limit of large outer 
scale it disappears from the final expressions. Since the 
power spectral density only depends on the magnitude 
of k, one can set k z = as discussed in Appendix [Bj so 
that k is replaced by K = | K| . 



A. Free-space propagation term 

First we consider the integral in Eq. (|3Tj) . which rep- 
resents the free-space propagation. After evaluating the 
integral and removing the superfluous mixed terms con- 
taining a p times a q, we obtain a generating function for 
the ^-integral, 



Sg(z) 



exp 



PrnPn 



+ w m )(l + w n ) 
8(1 - w m w n ) 3 
x [2(1 - w m w n ) + 



QmQrt 



2(1 - w m w n ) 



(37) 



where p m , p n , q m , q n , w m , and w n are the generat- 
ing function parameters, associated with the to- and n- 
indices. Since the p's (or g's) always appear in products 
for the to- and n-indices, the result implies an orthogo- 
nality condition for the azimuthal indices. The same is 
not true for the radial indices — one finds that the result 
is non-zero when the radial indices differ by either or 
1. When the difference is the final result for explicit 
modes can be expressed as 



i(l + \l\ + 2r) 
2z R 



(38) 
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where the azimuthal indices are indicated by l m — l n — I 
and the radial indices are indicated by r m = r n = r. 
When the difference between the radial indices is equal 
to f the result for explicit modes is given by 

_ ,< 1 + Vl + rpi + r) »i 

where the azimuthal indices are again indicated by l m = 
l n = I, but now the radial indices are indicated by r = 
(r m +r n - l)/2. 



B. Divergent dissipative term 

Next we consider the integral for the dissipative term 
given in Eq. (|3"2"]). Substituting Eq. into Eq. 
one obtains, 



0.1244C 2 A 



-2 -5/3 



(40) 



Hence, Lt diverges for large outer scale (k — > 0). How- 
ever, we find that these terms are canceled by similar 
terms coming from L m>n>u>v (z). 



C. Modal correlation functions 

The integral in Eq. (l34l) represents the correlation be- 
tween the momentum space wave functions of the OAM 
modes. To evaluate this integral we express the gen- 
crating function of Eq. (135[) in polar momentum space 
coordinates, so that k x + ik y = Kexp(i(f>). The result 
of the integration, left as a generating function, is then 
given by 



2(1 - w m w n ) 



exp 



PrnPn T" QmQn 
2(1 - W m W n ) 



iK exp(i<p)(p m Q + q n ( m )y 

2(1 - w m w n ) 
iK exp(-icj))(p n ( m + q m Q)v 
2(1 - w m w„) 

K 2 C m Qv 2 



2(1 - w m w n ) 



(41) 



where Cx — z r — iz — w x (zr + iz) and rj = X/lo . If we 
evaluate the azimuthal indices explicitly while leaving the 
radial indices implicit in terms of the parameter w m and 
w n , one can express the correlation function as 



W rG (K,^z) 



exp(-X) exp[i(m - n)(f)]E^ 



(1- 


- W m 


Wn) 




1/2 




[ r n ! " 


r 1 


_(\n\ + r n )l_ 




_(\m\ + r m )\ 



X 

M(m,n) 

x E 

s=0 



\m\l\n\l{-X)- s 
(|m| - s)!(|n| - s)!s! 



-,1/2 



(42) 



where m(— l m ) and n(= l n ) are the azimuthal indices; 
and r„ are their associated radial indices; and 



M{m,n) = ~(\m\ + \n\-\m-n\) (43) 



X = 
E m — 

E n = 



2(1 - w m w n ) 

iKCmV 
a/2(1 - w m w n ) 



(44) 
(45) 
(46) 



\/2(l - w m w n ) 

Note that W r a still represents a generating function with 
respect to the radial indices. 

D. General dissipative term 

The two-dimensional integration in Eq. (|33[) can be 
separated into a radial and angular integral in momen- 
tum space polar coordinates. Since &i(K) only depends 
on the radial coordinate (K = |K|) the integral over cj> 
only involves the (^-dependencies of the W-"s, as expressed 
in Eq. (|4"2"|. The combined ^-dependencies of the two 
W's is given by exp[i(m — u — n + v)<j>], where m, n, 
u and v are the azimuthal indices of all the modes in- 
volved. The result of the angular integral is zero unless 
in — u — Tt-\-v = 0, in which case the result is the product 
of the two W's without the ^-dependent exponentials, 
times 2n. As a result many of the elements in L mnuv 
vanish. 

The result of the remaining if-integral is too compli- 
cated to express as a single closed form expression. How- 
ever, one can consider the result on a term- by term basis. 
These terms all have the form 

Aexp{-BK 2 )K 2m 



f m (K) = 



(47) 



(iY 2 + K 2)H/6 ' 

where m is a non-negative integer (not to be confused 
with the combined OAM indices used before), A contains 
all the multiplicative parameters from Eqs. (|42|) and (|36[) . 
and B is a parameter composed of the parameters in the 
exponent of Eq. (|4"2")l . as given in Eq. (l4*4"l) . 

If m = the integral over f m (K) diverges as kq — >• 0. 
In the limit of small kq the leading terms are, 



f (K)K dK Lt 



rV2 



c 2 4 /3 (i + t^ 



6r(2/3) 



A 2 



(48) 

where t = z/zr and Lt is given by Eq. (|4T)|) . It turns out 
that the Lt term in Eq. exactly cancels all the Lt 
terms that appear inside the L mnpq term in Eq. (|30)) as 
a result of Eq. (|48|) . 

The integrals of f m (K) with m > all give finite re- 
sults independent of kq in the limit where kq — > 0, and 
have the form, 



fm(K)K dK « G 



A 2 



(49) 
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where G m is a numerical constant that only depends on 
m. 

VII. QUBIT EXAMPLE 
A. Solving the IPE 

Here we work through a simple example where only 
modes of the lowest radial index (r = 0) and with az- 
imuthal indices of the same magnitude I = ±q are con- 
sidered. We consider three cases where q = 1,2, 3, respec- 
tively. The small number of modes (only two per case) 
imply a severe truncation. The trace of the truncated 
density matrix is not in general equal to 1. The trun- 
cated density matrix can be normalized to calculate the 
concurrence [Til [l5j as a measure of the entanglement. 
On the other hand, the reduced trace gives an indication 
of the loss of information to the higher order modes. 

For r = one sets wi = W2 = in all the generating 
functions. Then all the non-dissipative terms (those that 
contain S) in Eq. (|3T)]) are equal and cancel each other, 
so that only the dissipative terms remain. 

After evaluating the integrals for L mi n iUiV (z) one finds 
that, in the limit of large outer scale, the only nonzero 
elements are 



q,q,q,q( z ) ~ L q q q q {z) 
9i9>9>9(' z ) = -Lq,q,q,q( z ) 



A q h{z) (50) 



and 



J q,q,q,q\ 



(z) = Lq t q !qtq (z) = B q h{z), 



(51) 



where q = —q: the quantities A q and B q are positive 
constants that only depend on q (A\ — 0.03976, B\ = 
0.0007675, A 2 = 0.05588,^2 = 0.0001213, A 3 = 
0.07110, B 3 = 0.00004444); and h(z) is the same func- 
tion for all the elements and contains all the dimension 
parameters, 



h(z) = 



i*\ (A. 
\7rw0 



2 \ 5/6 



(52) 



where zr is the Rayleigh range (ttluq/X), uj$ is the radius 
of the beam waist and A is the wavelength. 

The elements in Eq. (1501) . which contain the diverg- 
ing Lt from Eq. (|40[) . are the diagonal elements of 
L m ,n,u,v(z) when treated as a 4 x 4 matrix. One can 
therefore view Lt in Eq. (I50[) as being multiplied by an 
identity matrix. The Lt term in Eq. (|30p also represents 
an identity matrix, but with the opposite sign. This im- 
plies that the Lt term in Eq. (|30|) exactly cancel the 
Lt elements in Eq. ([50]) for both L m>n ,u,v(z) terms in 
Eq. (|3"0")l . leaving the final expression without Lt- As a 
result the outer scale drops out of the final expression. 

Assuming that the initial state of the density ma- 
trix is the singlet Bell state in the OAM basis — 



\q)\q))/ y/2, one obtains the following solution of the trun- 
cated density matrix 



T 



1 - R 2 

1 + R 2 -2R 

-2R 1 + R 2 

1 -R 2 



(53) 

where mp (nq) denote the row (column) indices, and 
where 



T = exp 
R = exp 



-2{A q - B q ) f h(z') dz 
Jo 

-2B q [ h(z') dz' 
Jo 



(54) 
(55) 



The eigenvalues of the truncated density matrix in 
Eq. ([53]) are T(l + R) 2 /4, T(l-R) 2 /4, T(l-i? 2 )/4 and 
T(l — i? 2 )/4, which are all positive. The trace of the trun- 
cated density matrix is given by T, which is a decaying 
function, since A q > B q . The trace indicates how much 
of the information is lost due to coupling to higher order 
modes that are not represented in the density matrix. 

The amount of entanglement for a bipartite qubit sys- 
tem is quantified by the concurrence [lj, [l5[ . Normaliz- 
ing the density matrix by setting T = 1, one can show 
that the concurrence for this case is given by 



C = ^(2R + R 2 - 1). 



B. Fast decay limit 



(56) 



One can evaluate the integral of h(z) in Eq. (|52|) . using 



[i + ttf* dt = wm(iv^-3) r /2 



x(l + i 2 )- Q-_\\% (it) 

64tt \3 
t + 0(t 3 ) 



(57) 



where t = zj zr and Q™ is an associated Legendre func- 
tion of the second kind [16| . For small t the result of the 
integral in Eq. (|57|) is approximately equal to t. In Fig. [2] 
the result in Eq. (|57|) is compared with the line given 
by t. One can see that for t < 1/3 the result of Eq. (|57| 
could be fairly well approximated by t. What this implies 
is that for propagation distances much shorter than the 
Rayleigh range the integral of h(z) can be approximated 
by, 



[ Z h(z') dz'« 0.5928 ( — 
Jo \ r o 



5/3 



(58) 



where r$ = 0.185(A 2 /C 2 /,z) 3 / 5 , which is the Fried param- 
eter. Thus all the dimension parameters are combined 
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into uja/rQ. If the entanglement completely decays over 
distances much shorter than the Rayleigh range — a sit- 
uation, which we call the fast decay limit — one can use 
Eq. (1551) and thus express the entire behavior simply as 
a function of u>o/fo- For larger values of t the behavior 
becomes more complicated as shown in Eq. ([57]) . In that 
case the behavior does not only depend on ujq /ro but also 
on the equivalent of ujo/tq with z replaced by zr. 




0.2 0.4 0.6 0.8 1 

Normalized propagation distance 



FIG. 2: Comparison of the t and the integral of h(z). 



C. Comparison 

Previous analyses of the decoherence of OAM entan- 
glement due to atmospheric scintillation [9l-[Tl| obtained 
results that only depend on ojq/tq. So to compare our re- 
sults with their results we need to consider our results in 
the fast decay limit. This is done by substituting Eq. (|55|) 
into R and T in Eqs. (|54|) and (|55|) . and then substitute 
R into Eq. ([56| . We plot the resulting trace T and the re- 
sulting concurrence C as functions of too/ro for q = 1, 2, 3 
in Figs.[Ha) and[3{h), respectively. 



(a) (b) 




FIG. 3: Plots of (a) the trace of the truncated density matrix 
and (b) the concurrence for a biphoton, initially in the singlet 
Bell-state, in terms of two OAM states with / = ±g, for q — 
1, 2, 3, as a function of uo/ro. 

From the curves for the trace in Fig. EHa) one can see 
that modes with higher OAM are scattered more rapidly 
into other modes than those with lower OAM. On the 
other hand, from Fig. [3th) we see that modes with higher 
OAM retain their entanglement for longer distances than 
those with lower OAM. These conclusions agree qualita- 
tively with previous work Q, however, the scattering into 



other modes occurs at a scale where ojo/tq ^> 1 and the 
entanglement lasts for at least two orders of magnitude 
longer than was found in Ref. Here the slowness of 
the decay in the concurrence is a result of the smallness 
of the values of the BqS. (While the A q 's come from the 
autocorrelation functions of the OAM modes, the B q 's 
come from the cross-correlations of modes with different 
azimuthal indices, which give much smaller overlaps with 
the power spectral density.) From these results it appears 
that the effect of scattering and the implied loss of pho- 
tons in the desired OAM modes may turn out to be a 
more significant challenge for free-space quantum com- 
munication than the decoherence of OAM entanglement. 
On the other hand, the effects of the severe trunction 
may imply that these results are but a poor reflection of 
what really would happen in this scenario. 



VIII. DISCUSSION 

The loss of entanglement due to decoherence is a chal- 
lenge that confronts the development of quantum com- 
munication systems. It is necessary to be able to predict 
the propagation scale over which this decoherence takes 
place before a successful free-space quantum communi- 
cation can be designed. Previous attempts to make such 
predictions [t| were based on certain assumptions Q , the 
effects of which were perhaps not completely clear. Intu- 
itively, the assumption that one can represent the atmo- 
sphere by a single phase screen sounds reasonable, and so 
does the assumption that one can represent the strength 
of the turbulence by the Fried parameter. However, care- 
fully considering the effect of turbulence, one realizes that 
some physical effects of the scintillation process is lost as 
a result of these assumptions. 

Turbulence is a cascaded process. The random index 
fluctuations cause a phase modulation of the traversing 
optical beam. Directly after an initial random phase 
modulation, the amplitude of the optical beam is un- 
affected. However, during subsequent propagation the 
phase distortion is partially transferred into an amplitude 
distortion. This mixture of phase and amplitude dis- 
tortion receives further phase modulations as the beam 
passes through the random medium. The random phase 
modulation and the propagation both play important 
roles in the process that produces the scintillated beam. 
Without the propagation part of the process the phase 
modulation will never be converted into amplitude scin- 
tillation. Since OAM modes have particular phase and 
amplitude characteristics, both the phase and the ampli- 
tude are important to give the correct scattering coeffi- 
cients. 

The formulation that is presented in this paper takes 
into account both the random phase modulation and the 
propagation and therefore gives the required effect on the 
phase and the amplitude of the modes. 

The effect of having a single phase screen where the 
turbulence is completely characterized by the Fried pa- 
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rameter is shown by the fast decay limit, discussed in 
Section IVIIBI Only in this limit can the effect of the 
turbulence be completely described by the Fried param- 
eter. Beyond this limit the scale parameters in the prob- 
lem starts to contribute in a way that cannot be com- 
bined into the Fried parameter. In other words, different 
combinations of these scale parameters can give different 
predictions even when the Fried parameter remains fixed. 
The implication is that, when one models the quantum 
scintillation process by a single phase screen in terms of 
the Fried parameter one tacitly assumes the fast decay 
limit and it is debatable whether the fast decay limit is 
valid or even useful for a free-space quantum communi- 
cation system. 



IX. SUMMARY 

We use the transformation of the momentum space 
wave functions of OAM modes after an infinitesimal 
propagation through a turbulent atmosphere to derive 
an IPE for the decoherence of OAM entangled biphoton 
states. The resulting set of first order differential equa- 
tions is used to study the evolution of a severely trun- 
cated density matrix where the initial state of a biphoton 
is an entangled qubit in the OAM basis. The results are 
compared with previous results in the literature. 



Appendix A: Generating function for OAM modes 

The Laguerre-Gaussian (LG) modes, which are solu- 
tions of the paraxial wave equation, are given in terms of 
normalized coordinates by 



M r L f 



. .Au±ivf\{i + ity 

,„ ( 2{u 2 +v 2 ) 
1 1 1 • t 2 



it -I 



(Al) 



where the normalized coordinates are given by u — x /uj$ , 
V = y/uJo and t — z / 'zr — zX/ttlo 2 , in terms of the initial 
radius of the mode profile ujq and the Rayleigh range zr\ 
r is the radial index (a non- negative integer); I is the 
azimuthal index (a signed integer); the ± sign is given 
by the sign of I; Af is a normalization constant given by 



H2IH+ 1 
<r+\l\) 



1/2 



(A2) 



and Li} represents the associate Laguerre polynomials, 
which can be obtained from the generating function, 



g w {x,w) 



1 



(i-wy+w 

by computing its r-th derivative 



exp 



—xw 
1 — w 



r rl dw r 



(A3) 



(A4) 



The azimuthal index / of the LG modes represent the 
amount of OAM that each photon in such an optical 
modes carries. Therefore, the LG modes act as an OAM 
basis in quantum optics. The LG mode functions are 
the position space wave functions of the OAM states and 
their Fourier transforms are the corresponding momen- 
tum space wave functions. 

In the paraxial limit the LG modes are treated as two- 
dimensional functions of the transverse coordinates u and 
v, and these two-dimensional functions change as a func- 
tion of the normalized propagation distance t. 

One can use a generating function for the OAM ba- 
sis functions (LG modes) to evaluate the integrals that 
contain OAM momentum space wave functions. Such in- 
tegrals only have to be solved once and afterward the 
solutions for particular OAM modes are generated using 
derivatives. Since it is always easier to compute deriva- 
tives than to solve integrals this represent a reduction in 
the computation effort that needs to be performed. 

The generating function for the LG modes in normal- 
ized coordinates is given by 



G = 



\ " 



1 



— T r ' 
m \ r, 



2(u 2 + v 2 ) 



n,m— 

[(u + iv)p + (u — iv)q] 



w(l + it) 



I -it 



1 



Q(t,w 

(l + w)(u 



(i - ity+ m 

(u + iv)p (u — iv)q 
Q(t,w) 



■ exp 



i „,2 



Sl(t,w) 



(A5) 



where f2(i, w) — 1 — w — it — iwt. The parameters p, q 
and w are used to generate particular LG modes in the 
following way, 



M r L f(u,v,t) 



for I > 

w 7 p 7 q— 

\dZ,G] n for ; = 

for / < 0, 



hd r J?G 



w,p,q—0 



(A6) 

where r and I represent the radial and azimuthal indices, 
respectively, and f\f is the normalization constant given 
in Eq. (|A"2|) . 

The integrals in Eqs. (|3T1) and ([34]) contain the Fourier 
transform of the LG modes. For this purpose we need 
the Fourier transform of the generating function, which 
is given by 



T{G} 



1 



exp 



iir(a + ib)p iir(a — ib)q 



1 



1 



7T 2 (a 2 + b 2 )n(t lW ) 



1 



(A7) 



where a and b are normalized spatial frequency compo- 
nents that are related to k x and k y through 

2wa 2irb 

CJ UJQ 



(A8) 
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The Fourier transforms of particular LG modes are ob- 
tained using the same procedure given in Eq. (|A6|) . 



autocorrelation function of the random function is given 
by 



Appendix B: Ensemble average 

As mentioned in Sec. IIII1 the refractive index fluctua- 
tions produced by a turbulent atmosphere are small com- 
pared to the average refractive index of air, Sn = n — 1 <C 
1, which leads to the fact that one can separate the prop- 
agation through a turbulent atmosphere into two parts: 
free-space propagation and the random phase modula- 
tion. The random phase function for the latter step is 
obtained by integrating the refractive index fluctuations 
through a thin sheet of atmosphere along the propagation 
direction, 

/•Zo+Sz/2 

0(x,y) — k / Sn(x, y, z)dz 

Jza-Sz/2 

fa k 8z 5n{x,y, z$), (Bl) 

where, in the last line we took the limit Sz —> 0. Re- 
placing the refractive index fluctuation with its Fourier 
expansion, we obtain 



(at, y, z ) = kSz j exp[—i(k x x + k y y + k z z )] 
c\ 3 k 

xNjk) 



(2tt) 3 ' 



where N n (k) is the three-dimensional spatial spectrum 
of index fluctuations. One can now define a two- 
dimensional spectrum for the accumulated index fluctu- 
ations over a thin sheet of atmosphere as follows, 



N n {K,z) 



f dk 

/ exp(-ik z z)N n (k) — i, (B3) 

J 2lT 



where A„(K, z) is the two-dimensional spectrum, which 
depends on the z position of the thin sheet. The three- 
dimensional spectrum of the refractive index fluctuations 
can be expressed in terms of its three-dimensional power 
spectral density, which follows from the correlation func- 
tion of the index fluctuations and which represents the 
model for the turbulence, 



iV n (k)=x„(k) 



$ (k) 



A 3 fc 



1/2 



(B4) 



where x(k) is a normally distributed random complex 
spectral function and Afc is its spatial coherence length 
in frequency domain. The latter is inversely proportional 
to the outer scale of the turbulence. Since the refractive 
index fluctuation Sn is an asymmetric real-valued func- 
tion, we have that X*(k) = £(— k). Furthermore, the 



(X(ki)x*(k 2 )) - (2^A fc r <5 3 (k! - k 2 ). (B5) 
In Eq. (f2"Tj) we find ensemble averages inside double 
^-integrals. Substituting Eqs. (IB3|) and (IB4[) into these 
ensemble averages, using Eq. (|F35[) to evaluate the ensem- 
ble average, one obtains 



(N n (K 1 ,z 2 )N*{K 2 ,z 1 )) dz 2 d Zl 



Zq J Zq 
2 i 



(2t:) z S(K 1 



K 2 



*o(ki) 

dk z 



x exp [ik z (zi - z 2 )] dz 2 dz 1 — — 

Z7T 



(B6) 



Setting z = zq + dz, we evaluate the two z-integrals 

pZQ-\-dz pzi 

/ / exp [ik z (zi — z 2 )\ dz 2 dzi 

J Zq J Zq 



1 



cos(k z dz) ,sm(k z dz) — k z dz 



1.2 

A. y 



(B7) 



The power spectral density $g(ki) is always even in k z . 
Therefore, the imaginary part of Eq. (|B7j) . being odd in 



(B2) 

have 



k z , does not contribute to the final expression. So we 



f f \N n (K u z 2 )N*(K 2 , Zl )) dz 2 dzi 

J Z J Z 

(2tt) 2 S(K 1 - K 2 ) J * (ki) 



1 — cos(k z dz) 



dk z 
2tt 



(B8) 



Due to the fact that the refractive index variations are 
very small, the light that propagates through the turbu- 
lent atmosphere remains unchanged over distances much 
longer than the correlation distance of the turbulent 
medium. One can therefore assume that dz is much 
larger than this correlation distance. As a result the 
function inside the square-brackets in Eq. (|B8|) acts like 
a Dirac delta function, so that one can substitute k z — 
in $o and pull it out of the fc z -integral. The integral can 
then be evaluated to give 

f 1 f (N n (K u z 2 )N*(K 2 ,z)) dz 2 dz 

J Zq J Zq 

= 2 7 r 2 dzS(K 1 -K 2 )*i(Ki), (B9) 



where we defined <i>i(Ki) = $o(Ki,0). 
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